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We investigate theoretically the simultaneous tunneling of two electrons from a superconductor 
into a normal metal at low temperatures and voltages. Such an emission process is shown to be 
equivalent to the Andreev reflection of an incident hole. We obtain a local tunneling Hamiltonian 
that permits to investigate transport through interfaces of arbitrary geometry and potential barrier 
shapes. We prove that the bilinear momentum dependence of the low-energy tunneling matrix 
element translates into a real space Hamiltonian involving the normal derivatives of the electron 
fields in each electrode. The angular distribution of the electron current as it is emitted into the 
normal metal is analyzed for various experimental setups. We show that, in a full three-dimensional 
problem, the neglect of the momentum dependence of tunneling causes a violation of unitarity and 
leads to the wrong thermodynamic (broad interface) limit. More importantly for current research 
on quantum information devices, in the case of an interface made of two narrow tunneling contacts 
separated by a distance r, the assumption of momentum-independent hopping yields a nonlocally 
entangled electron current that decays with a prefactor proportional to r~ 2 instead of the correct 
r-\ 
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I. INTRODUCTION 



The electric current through a biased normal- 
superconductor (NS) interface has for long been the ob- 
ject of extensive theoretical and experimental attention 
HI, 0, Lj- lH ■ Recently, new interest in this classic problem 
has been spurred by the possibility of using conventional 
superconductors as a natural source of entangled electron 
pairs that may be injected into a normal or ferromagnetic 

metal [1 EH Ell El El El El El El El El El and 

eventually used for quantum communication purposes. 
Clearly, the efficient and controlled emission of electron 
singlet pairs into normal metals or semiconductor nanos- 
tructures requires a deeper understanding of the underly- 
ing transport problem than has so far been necessary. In 
particular, it is of interest to investigate how the entan- 
gled electron current depends on various parameters such 
as the shape and size of the NS interface as well as the po- 
tential barrier profile experienced by the tunneling elec- 
trons. A preliminary focus on tunneling interfaces seems 
adequate, both because such interfaces are amenable to 
a simpler theoretical study and because the low electric 
currents which they typically involve will facilitate the 
control of individual electron pairs. 

In the light of this new motivation, which shifts the 
attention onto the fate of the emitted electron pairs, it 
seems that the picture of Andreev reflection, which so 
far has provided an efficient book-keeping procedure, has 
reached one of its possible limits. When dealing with fi- 
nite size tunneling contacts, the Hamiltonian approach 
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is more convenient than the calculation of the scatter- 
ing wave functions, since it does not require to solve the 
diffraction problem to find the conductance. Moreover, 
it is hard to see how problems such as the loss of non- 
local spin correlations among distant electrons emitted 
from a common superconducting source can be analyzed 
in terms of Andreev reflected holes in a way that is both 
practical and respectful to causality. While an Andreev 
description may still be practical in situations involving 
multiple electron-hole conversion, the fate of the quasi- 
particles in the outgoing scattering channels will have to 
be investigated in terms of a two-electron (or two-hole) 
picture if one is interested in studying nonlocal correla- 
tions in real time. 

Recently, several authors 0H H EE O El, EE E3, El 

have addressed the emission of electron pairs through two 
distant contacts in a language which explicitly deals with 
electrons above the normal Fermi level. Here we investi- 
gate the emission of electron pairs from a superconductor 
into a normal metal through tunneling interfaces of differ- 
ent geometrical shapes and potential barriers. With this 
goal in mind, we devote section II to rigorously estab- 
lish the equivalence between the pictures of two-electron 
emission and Andreev reflection of an incident hole. We 
argue extensively that each picture reflects a different 
choice of chemical potential for the normal metal, a point 
also noted in Ref. 19] . After a precise formulation of 
the problem in section III, we derive a real space tun- 
neling Hamiltonian in section IV that accounts for the 
fact that electrons with different perpendicular energy 
are transmitted with different probability through the 
interface. In section V, we study the structure of the 
perturbative calculation that, for vanishing temperatures 
and voltages, will yield the electron current to lowest or- 
der in the tunneling Hamiltonian. Section VI concerns 
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FIG. 1: Hole Andreev reflection vs two-electron emission: (a) 
When fj,s is used as the reference chemical potential both in 
S and N, a typical scattering process at the NS tunneling in- 
terface for eV = fj,s — hn > is that of incoming holes with 
energies between [is and /is + eV that are most often nor- 
mally reflected and only rarely Andreev reflected, (b) If, al- 
ternatively, /in is used to define quasiparticles in N, the many 
incoming or outgoing holes are viewed as a vacuum of quasi- 
particles. The outgoing electron generated in a rare Andreev 
reflection event appears now as a spontaneously emitted elec- 
tron above us- Such an event causes an outgoing hole state to 
be empty. This is now perceived as the emission of a second 
electron with energy between /ijv and us- Tracking the spin 
and the momentum component parallel to the interface leads 
to the picture of two electrons emitted in a spin singlet state 
with opposite parallel momenta. 



itself with the angular dependence of the current through 
a broad NS interface, providing the connection with cal- 
culations based on the standard quasiparticle scattering 
picture pfl l2l| . In section VII we investigate the tun- 
neling current through a circular NS interface of arbi- 
trary radius, paying attention not only to the total value 
of the current but also to its angular distribution and 
to the underlying two-electron angular correlations. We 
also investigate how the thermodynamic limit is achieved 
for broad interfaces. Section VIII deals with the elec- 
tron current through an interface made of two distant 
small holes, focussing on the distance dependence of the 
contribution stemming from nonlocally entangled elec- 
trons leaving through different holes. In section IX we 
investigate the commonly used energy-independent hop- 
ping model and prove that it violates unitarity, leads to 
a divergent thermodynamic limit, and yields the wrong 
distance dependence for the current contribution coming 
from nonlocally entangled electrons. A concluding sum- 
mary is provided in section X. 



II. TWO-ELECTRON EMISSION VS. ANDREEV 
REFLECTION 

In a biased normal-superconductor tunneling interface 
in which e.g. the superconductor chemical potential is 
the greatest, one expects current to be dominated by 
the injection of electron pairs from the superconductor 
into the normal metal if the voltage difference V and the 
temperature T are sufficiently low, single-electron tun- 



neling being forbidden by the energy required to break a 
Cooper pair. Specifically, one expects two-electron tun- 
neling to dominate if fcgT, eV <C A, where A is the 
zero-temperature superconductor gap. Simple and un- 
questionable as this picture is, it is not clear how it can 
be quantitatively described within the popular Bogoli- 
ubov - de Gennes (BdG) quasiparticle scattering picture 
0iiil. While it leaves the BCS state unchanged, 
the emission of two electrons into the normal metal in- 
volves the creation of two quasiparticles, something that 
is not possible within the standard BdG formalism, where 
the quasiparticle number is a good quantum number and 
the quasiparticle scattering matrix is thus unitary. The 
conservation of quasiparticle current is a consequence of 
the implicit assumption contained in the conventional 
BdG scheme that the reference chemical potential used 
to identify quasiparticles in the normal metal is the su- 
perconductor chemical potential /is. However, as shown 
below, one does not need to be constrained by such a 
choice. 

In the mean field description of inhomogeneous super- 
conductivity provided by the BdG formalism, the Hamil- 
tonian is given by 



H = E 



^ 1 Tno-Tno-) 

ii n 



(1) 



where E is the condensate energy and 7^ creates a 
quasiparticle of energy e„, spin quantum number a and 
wave function [u n (r),v n (r)] satisfying the BdG equations 
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where Hq = — ft 2 V 2 /2m + U is the one-electron Hamil- 
tonian. In the standard convention one adopts solutions 
such that s n — /is > 0. However, a fundamental property 
of the BdG equations |22l l24l |25| is that, for every quasi- 
particle n, a of energy e n > MSj there exists another solu- 
tion n', a' with spin a 1 = —a, energy e n > = —e n + 2/zg < 
Us and wavefunction ( ) = (-«*,«*). These two 

solutions are not independent, since creating quasiparti- 
cle n,a is equivalent to destroying quasiparticle n',—a 
[23. More specifically, 7^ = 7„/ T , and 7„ T = -7,^. 

In the case of a normal metal, where quasiparticles are 
pure electron or pure holes, the above property implies 
that creation of a quasiparticle of energy e n > Ms and 
wave function (0, v n ) (i.e. a pure hole) corresponds to the 
destruction of a quasiparticle of energy e' n — 2/i^ — e„ < 
jj,s and wave function (— d*,0) (i.e. a pure electron). If 
v n (r) ~ exp(ik/ l -r), the existence of a hole of momentum 
k/j, with kh < hp, and energy e > fj-s corresponds to 
the absence of an electron in the state of wave function 
v*(r) ~ cxp(-ik,j • r) with energy e' = 2pt s - £ < Us- 

In a biased NS tunneling structure, the normal metal 
has a different chemical potential /ijy = /J-s — eV. With- 
out loss of generality, we may assume /in < Ms- If we 
release ourselves from the standard BdG constraint of 
using [xg as the reference chemical potential even on the 
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normal side, a clearer picture is likely to emerge. We 
may decide that, in the energy range /Ltjv < s' < /is, we 
switch to the opposited convention for the identification 
and labelling quasiparticles. In other words, we decide to 
use Hn as the reference chemical potential. Translated 
to the example of the previous paragraph, we pass to 
view the occupation of the hole- type quasiparticle state of 
wave function (0, v n ) and energy e > (As as the emptiness 
of the electron-type quasiparticle state of wave function 
(— v*,0) and energy e' < fig. Conversely, the absence of 
quasiparticles in (0, v n ) is now viewed as the occupation 
of (— «*, 0), i.e. as the existence of an electron with wave 
function — w*(r) and energy e' between /zjy and (is- 

The consequences that this change of paradigm has 
on the way we view transport through an NS interface 
can be more clearly appreciated in Fig. ^ In the stan- 
dard BdG picture represented in Fig. [TJa), with (is as 
the reference chemical potential, the "in" state is that 
of many holes impinging on the NS interface from the N 
side, with energies between (is and (As+eV. Since ours is 
a tunneling structure, normal reflection is the dominant 
scattering channel and only one hole is Andreev reflected 
as an electron (quasiparticle transmission is precluded at 
sufficiently low temperatures and voltages). Thus the 
"out" state is that of many holes and only one electron 
moving away from the surface, all with energies also be- 
tween (is and (As + eV, since quasiparticle scattering is 
elastic. Given the unitary character of quasiparticle scat- 
tering in the BdG formalism, the existence of an outgoing 
electron requires the outgoing hole quasiparticle state at 
the same energy to be empty, due to the incoming hole 
that failed to be normally reflected. The absence of such 
an outgoing hole is clearly shown in Fig. ^a). 

If we now shift to (An as the reference chemical po- 
tential, the picture is somewhat different. The many 
impinging holes on the surfaces are now viewed as the 
absence of quasiparticles, i.e. the "in" state is the vac- 
uum of quasiparticles. The one electron that emerged 
from a rare Andreev process continues to be viewed as 
an occupied electron state, shown above (is in Fig. ^b). 
The many outgoing holes of the BdG picture are again 
viewed as an absence of quasiparticles. The second out- 
going electron that is needed to complete the picture of 
two-electron emission corresponds to the empty outgo- 
ing hole state of the BdG picture which originates from 
the hole that failed to be normally (and thus specularly) 
reflected. It is shown in Fig. ^b) with energy between 
(In and (is- As is known from the theory of quasiparti- 
cle Andreev reflection, the outgoing electron of Fig. QJa) 
follows the reverse path of the incident hole (conjugate 
reflection). Therefore the two electrons in Fig. ^b) have 
momenta with opposite parallel (to the interface) compo- 
nents and the same perpendicular component, i.e. they 
leave the superconductor forming a V centered around 
the axis normal to the interface. Inclusion of the spin 
quantum number completes the picture of two electrons 
emitted into the normal metal in an entangled spin sin- 
glet state. 
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FIG. 2: Schematic lateral view of the NS tunneling structures 
studied in this paper: (a) circular interface of arbitrary radius 
R and (b) interface made of two small holes at a distance r 
from each other. The rest of the NS interface is assumed to 
be opaque. 



In summary, we have rigorously established the equiva- 
lence between the pictures of Andreev reflection and two- 
electron emission, noting that they emerge from different 
choices of the chemical potential to which quasiparticles 
in the normal metal are referred, (is in the standard BdG 
picture, and (An in the scenario which contemplates the 
spontaneous emission of two electrons. For simplicity, 
and because it better fits our present need, we have fo- 
cussed on the case of a tunneling structure. However, 
the essence of our argument is of general validity. Here 
we just note that, in the opposite case of a transmissive 
NS interface |3|, |26|, the same argument applies if, ex- 
changing roles, Andreev reflection passes to be the rule 
while normal reflection becomes the exception. In that 
case, charge accumulation and its accompanying poten- 
tial drop, which are generated by normal reflection |27j . 
will be essentially nonexistent. 

Upon completion of this work, we have learned that 
the need to change the normal metal vacuum to describe 
hole Andreev reflection as electron emission has also been 
noted in Ref. Gjj|. 



III. FORMULATION OF THE PROBLEM 

As has been said, an extensive body of literature has 
been written on the various aspects of electron transport 
through a normal-superconductor interface El 13, JUT 



[m El M El m M IS M El ElHIIli 

I3d I3f| . Generally those works have focussed on the case 
of broad interfaces or point contacts [2! . Our goal 
here is to analyze the current of spin entangled Cooper 
pairs from a BCS bulk superconductor into a bulk normal 
metal through an arbitrarily shaped insulating junction 
in the tunnel limit. Apart from the desire to explore 
novel types of NS structures, we are also motivated by 
the need to investigate in depth the two-electron emis- 
sion picture, which is likely to be useful in the design of 
quantum communication devices. We wish to consider 
explicitly geometries of the sort depicted in Fig. [3 i.e. 
a 2D planar interface of arbitrary radius R, presented 
in Fig. Ufa), and two small orifices separated by a dis- 
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tance r, shown in Fig. Efb). It is assumed that, out- 
side the designed region, the interface is opaque to the 
flow of electrons. For simplicity, both the normal and 
the superconducting electrodes are taken to be ballistic. 
An advantage of the tunneling regime is that the prox- 
imity effect may be neglected, i.e. we assume that the 
gap function drops sharply at the NS interface and that 
self-consistency in the gap may be safely neglected pBj . 
Another benefit is that we deal at most with two chem- 
ical potentials, since the low scale of tunneling currents 
gua rantees that the normal metal is close to equilibrium 
|27| and that no phase slips develop within the supercon- 
ductor [3^1 . Inelastic processes at the interface will also 
be ignored |33j . 

We are interested in a conventional (s-wave) super- 
conductor because it may act as a natural source of 
spin-entangled electrons, since its electrons form Cooper 
pairs with singlet spin wave functions and may be in- 
jected into a normal metal. The superconductor, which 
is held at a chemical potential jj,s, is weakly coupled by 
a tunnel barrier to a normal metal which is held at fj.^. 
By applying a bias voltage V — (/is — /iAr)/e such that 
eV > 0, transport of entangled electrons occurs from 
the superconductor to the normal metal. We focus on 
the regime k B T < eV < A. Since 5 = A/E F ~ 10" 4 
in a conventional superconductor, rearrangement of the 
potential barrier due to the voltage bias can be also ne- 
glected. However, the effect of a finite, small S will often 
be tracked because pairing correlations (and thus non- 
local entanglement) decays on the scale of the coher- 
ence length £o> which is finite to the extent that A is 
nonzero. For convenience we assume that the supercon- 
ductor normal-state properties (m, k F , etc.) are the same 
as for an ordinary metal. 

We will use a tunneling Hamiltonian approach and ex- 




where zo lies inside the barrier, i.e. Jkq(r, zo) = (i/^-)^kq 
is the matrix element of the z component of the current 
density operator. Due to charge conservation, Jk q is in- 
dependent of the choice of point zq G [— w/2,w/2]. The 
unperturbed wave functions are of the form 

e ik||r 

Xk(r,z) = -j=-tp k „(z), (6) 

where the exact shape of ipk, (z) depends on the barrier 
height. Thus, 

Hereafter, the volume of each metal £Il,r is taken equal 
to = AL, A being the area of the interface and L 



plicitly consider the emission of two electrons from the 
superconductor, a viewpoint that will be mandatory in 
contexts where the late evolution of correlated electron 
pairs in the normal metal is to be investigated. 



IV. THREE-DIMENSIONAL TUNNELING 
HAMILTONIAN 

The Bardeen model for electron tunneling [34j assumes 
that a system made up of two bulk metals connected 
through an insulating oxide layer can be described by 
the Hamiltonian 

H = H L + H R + H T . (3) 

Here Hp and Hr are the many-body Hamiltonians for 
the decoupled (i.e. unperturbed) electrodes, the super- 
conductor being on the left and the normal metal on 
the right. The connection between both electrodes is de- 
scribed by the tunneling term Ht (see e.g. Ref. lii^l ;: 

ffT = ^T kq 4v+H.c. (4) 

kqtr 

Here is the creation operator in the normal metal 
of the single-particle state of orbital quantum number k 
and spin cr, whereas c q(T destroys state q, a in the super- 
conductor and Tk q is the matrix element connecting both 
states. We assume a perfect interface defined by a square 
barrier t/(r, z) = U <d(z + w/2)Q{w/2 - z)] (hereafter r 
refers to the in-plane coordinate). 

If Xq( r ; z ) are the left-side stationary waves for a po- 
tential step Ul{y,z) = Uo<d(z + w/2) and Xk(r, z) be- 
haves similarly for Un(r, z) = Uq®(w/2 — z), Bardeen 
PI showed 



(5) 

Z = Zq 

I 

the length of each semi-infinite metal. N(Q) is the 3D 
one-spin electronic density of states of the normal metal 
at the Fermi level: N(0) = k F /Ait 2 E F . We define the 
transparency of the barrier as 

r = 4<[?fe-v° w , (8) 
V Uq 

where pq = \/2mUo/h. In the particular case pow ^> 1 
and Ep <C Uq, t coincides with the probability amplitude 
that an electron with perpendicular energy E z = Ep 
traverses the barrier. L(k z ,q z ) in Eq. J7J) captures the 
dependence of the hopping energy on the z momentum 
component. Some authors take it as constant, but we 
shall argue in section IX that its k z , q z dependence is 
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crucial for a sound description of 3D transport problems. 

For a square barrier, we may define u = Uy/Ep, g z = 
Qz/kF, k-z = k z /kp, and write 

L(k z ,q z ) = k z q z a{n z ,g z ) exp{p w[l - a(K Z) g z )]} , (9) 

where 

a(x,y) = [b(x)+b(y)}/2, (10) 
b(x) = y / l-x 2 /u. (11) 

For high barriers (u — > oo) we have a(x) — > 1. Then, 

If we make Uq — > oo while keeping the electron transmis- 
sion probability finite, we are implicitly assuming that 
the barrier becomes arbitrarily thin (w — > 0), i.e. we are 
taking it to be of the form V{z) — H5(z), as popular- 
ized in Ref. On the other hand, since the height of 
the barrier is judged in relation to the perpendicular en- 
ergy E z < Ep < Uq, it is clear that, given Uq and w, Eq. 
(|12fl becomes correct for sufficiently small k z ,q z . In other 
words, Tk q behaves identically for u — ► oo or k z , q z — > 0. 
As a consequence, such bilinear dependence of Tk q for 
sufficiently small k z ,q z may be expected to hold for ar- 
bitrary barrier profiles within the tunneling regime. We 
note that Eq. (|12fl differs from the result obtained in Ref. 
[38| for the low energy hopping. 



Thus we adopt as a criterion of the validity of Bardeen's 
approximation that Eq. Ijlfijl holds, which from (|15|l . im- 
plies T(Ep) <C 1. This defines an upper energy cutoff E c 
in the various sums over electron states, which is the max- 
imum energy for which the approximation (|16l) is valid. 
For the square barrier, E c ~ Uo — h 2 /2mw 2 . 

For processes described by amplitudes which are first 
order in Ht, and as long as Uq is high enough compared 
to Ep to fulfill condition (|16|) for all relevant E z , all elec- 
tron momenta lie within the applicability of the tunnel 
limit and we may use the tunneling Hamiltonian safely. 
That is the case of the tunnel current through a NN in- 
terface or the quasiparticle tunnel current through a NS 
interface. 

The situation is different for transport through a NS 
interface, since it requires the coherent tunneling of two 
electrons. Then, the leading contribution to the tunnel- 
ing amplitude is quadratic in Ht and the final transmis- 
sion probability is sensitive to the existence of interme- 
diate virtual states where only one electron has tunneled 
and a quasiparticle above the gap has been created in 
the superconductor. Unlike the weighting factors of the 
initial and final states, which are controlled by the Fermi 
distribution function, the contribution of the virtual in- 
termediate states decays slowly with energy and the cut- 
off E c may be reached. In section V we show that there 
are two cases where the cutoff can be safely neglected, 
namely, the limit of high barrier («>1) and the limit 
of small gap (S <g; 1). 



A. Validity of the tunneling Hamiltonian model: 
momentum cutoff 

We wish to quantify the idea that a perturbative treat- 
ment of Bardeen's tunneling Hamiltonian is valid only 
when it involves matrix elements between weakly cou- 
pled states 0, |3(| . 

The transmission probability for a low energy electron 
incident from the left can be written 

T{E Z ) = WJJ< 1 , (13) 

where J q is the current density carried by the incoming 
component of the stationary wave q, and 

Wq^^E |7kq| 2 <5(£k - Eq) (14) 

k 

is the tunneling rate. Using Eqs. J7J) and ©, Bardeen's 
theory yields 

T(E Z ) = 16^ (l-^Ae- 2 ^, (15) 

where p z = ^/2m(Uo — E z )/h. On the other hand, an 
exact calculatio n 1371 recovers the tunneling result l|15(l 
for E z < Uo |39j, |4fj if we make the approximation 

sinh(p z uj) w cosh(p z w) w e PzW /2 . (16) 



B. Tunneling Hamiltonian in real space 

One of our main goals is to investigate transport 
through tunneling interfaces of arbitrary shape |4l| that 
are otherwise uniform. For that purpose we need a reli- 
able tunneling Hamiltonian expressed in real space. Our 
strategy will be to rewrite Eq. @ as an integral over the 
infinite interface and postulate that a similar Hamilto- 
nian, this time with the integral restricted to the desired 
region, applies to tunneling through the finite-size inter- 
face. The discontinuity between the weakly transpar- 
ent interface and the completely opaque region causes 
some additional scattering in the electronic wave func- 
tions that enter the exact matrix element. However, this 
effect should be negligible in the tunneling limit. In fact, 
we provide in Appendix A an independent derivation of 
the continuum results shown in this section which starts 
from a discrete tight-binding Hamiltonian. 

Thus, in we introduce the transformations 

c ta = J dr Jdz Xk(r, z) ^ N {r, z;a) (17) 

c qCT = J dr J dz x*(r, z) i() S (r,z;a) (18) 

where the wave functions x q and Xk are, respectively, 
solutions of Hl and Hr and are given in ©. ip\j(r, z; a) 
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and ips( r , z; a) are the field operator in the normal and Invoking Eq. (JTJ and the completeness of plane waves 
superconducting metals. in the x, y plane [which yields a term 5(r — r')], we obtain 



Ht = E ^^(o) j A dr j dz ' J dz L( - z > z ') ^( r ' z; CT ) z ' ; + Hx - ' ( 1Q ) 



where 



L(z,z') 



1 



E 



<p*(z')L(k x ,q g 



(20) 



Since the initial Hamiltonian (0} connects states which 
overlap in a finite region below and near the barrier, it 
is logical that the real space Hamiltonian l|19|l is non- 
local in the z-coordinate. An interesting limit is that of 
a high and (to keep transmission finite) thin barrier, i.e. 
the delta barrier limit. Then, the perpendicular wave 
functions can be precisely written 



Vkziz) — y/2/ Lsm(k z z), z > 



(21) 



and similarly for the left electrode. We introduce such 
wave functions in Eq. (|20|l and invoke the identity (here- 
after L — > oo ) 



— S~] k z sin(k z z) = -S'(z), 



k z >0 



(22) 



where the volume per orbital in fc z -space is tt/L. Then, 
to leading order in u^ 1 <C 1, Eq. (|19|l yields 



Hi 



E 



8^(0) J A 



dr 



d^ N (r,z;a) 



(hpsjr, z';a) 
dz' 



H.c. 



(23) 



If we replace the thermodynamic area A by a specific 
finite area, the real space Hamiltonian l(23|l can be used to 
describe tunneling through interfaces of arbitrary shape. 
As we have said, in Appendix A we provide an alter- 
native derivation which makes Eq. 1)23(1 appear as the 
natural continuum limit of the hopping Hamiltonian in 
a regularized tight-binding representation. We note that 
the tunneling Hamiltonian l(23|l may also be obtained if 
the r.h.s. of Eq. I(21|l is replaced by a plane wave repre- 
sentation. 

^From Eq. iffify we conclude that apparently reason- 
able choices of local tunneling Hamiltonian such as those 
oc J ifigipN lead to unphysical results in 3D. This point 
will be discussed in depth in section IX. 

To describe tunneling in real space, rather than 
starting from Hamiltonian (|19|1 . or its limiting version 
(|23|l . it is more convenient in practice to go back to 
Eq. |0 and make the replacement 5(kii — cjri ) — > 
(27r)- 2 / A dre i ( k !|- q n) r , with A finite. Then Eqs. O 
and H2Ufl may equivalently be written 

(24) 



If we make L(k z ,q z ) — k z q z , it is easy to prove that l|24|) 
becomes 



V. PERTURB ATIVE CALCULATION OF THE 
STATIONARY CURRENT 

Following Ref. we write the stationary electron 
current from the superconductor to the normal metal as 

Im = 2eJ2WfiPi, (25) 

if 

where Wfi is the transition rate at which two electrons 
tunnel from the superconductor into the normal metal, 
and pi is the stationary distribution accounting for a 
chemical potential difference eV between the two elec- 
trodes. We calculate the transition rate with a T-matrix 
approach |4^ |. 

Wf i = ^\{f\f( Ei )\i)\*8(e f -e i ). (26) 
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The T-matrix can be written as a power series in the 
tunnel Hamiltonian Ht, 

oo 

f (e) = H T + H T Y^[G (e)H T ] n (27) 
n=l 

where Go(e) = (e — i?o + «0 + ) _1 is the retarded Green 
function for the decoupled system. 

At zero temperature the initial state is \i) — |F) ® 
|BCS), where |F) is the filled Fermi sea ground state of 
the normal metal and |BCS) is the BCS ground state of 
the superconductor. The state \i) is the vacuum of quasi- 
particles if these are referred to fis in the superconductor 
and to /i at = /is — eV in the normal metal (see section 
II). In the final state 

i/) = ^(4 lT i,i-4i4T)i i ). ( 28 ) 

i.e. the superconductor remains unperturbed within the 
BCS description, since an entire Cooper pair has been re- 
moved, and two singlet-correlated electrons hover above 
the normal Fermi sea |43|. In the ensuing discussion we 
take e% = 2/is = 0. 

Since we wish to focus on the regime ksT <C eV <C A, 
single electron emission is forbidden due to energy con- 
servation, because it requires the breaking of a Cooper 
pair. Therefore, to leading order in Ht, we may approx- 
imate 

f (0) « T" = H T Go H t (29) 

and so we write 

(f\f (0)|i) = ((c k2i c klT - c k2T c kli ) T"> . (30) 

We insert a complete set of single-quasiparticle (vir- 
tual) states, i.e. 1 = £ kqCTCT , 7<U C L' I 7q<r, be- 
tween the two Ht in 1)29(1 and we use the fact that the 
resulting energy denominator |i0 + — £ k — Eq\ « |£q|, since 
£ k = h 2 k 2 /2m — /ijv — * when eV^ — * 0. To see this, one 
must note that the energy conservation implies £/ = 
therefore, £/ = £ kl + £ k2 — 2eV~ = 2[i s — 0. Thus, when 
eV — > 0, one may write £ kl w — £ k2 w . We also make 
use of MqVq = ii_ q u_q. Finally, we get 

= 2v / 2(c k2T cL T ^a4 li )^ (CqT ^ qi> r kiq r k2 ,_ q (31) 

where S q = [(fi 2 /2m) 2 (q 2 - fc 2 ,) 2 + A 2 ] 1 / 2 is the quasi- 
particle energy and F q = (c a TC_ a i ) is the condensation 
amplitude in the state q [22]. 

At zero temperature we have F q = A/2£q. Thus, in 
the summation of Eq. I|31l) , the contribution of high en- 
ergy virtual states is weighted by the Lorentzian Fq/Eq, 
of width A and centered around Ep. We already men- 
tioned in the previous section the need for a high-energy 




FIG. 3: Angular dependence of the normalized tunnel cur- 
rent ranging from ?7o /Ef — > oo for the outer (cos 5 6) curve 
to Uo/Ef = 1.1 for the inner one. Finite barriers have a 
width w = 5Af- Observe how the angular distribution focal- 
izes around the perpendicular direction as the barrier hight 
decreases. 

cutoff E c to prevent the inclusion of states for which the 
tunneling approximation is not valid. However, in the 
limit A/Ep — > 0, the Lorentzian distribution becomes 
a delta function and the results are independent of the 
cutoff, which can be safely taken to infinity. A similar 
situation is found in the limit Uq S> Ep, for which the 
sum in Eq. 131|) converges before reaching the energy E c 
above which Bardeen's approximation is no longer valid. 
In any of these two limits (A/Ep, Ep /Uq <§C 1), it is 
correct to take E c — > oo. 

VI. TOTAL CURRENT AND ANGULAR 
DISTRIBUTION THROUGH A BROAD 
INTERFACE 

The current through a NS junction is most easily cal- 
culated when the interface section is much bigger than 
A p . We shall refer to it as the broad interface or thermo- 
dynamic limit. Its detailed understanding is of interest 
for later reference in the investigation of finite size inter- 
faces. For ksT <C eV <C A, the set of equations in the 
previous section yields 

Iks=Iv-z- d6 sinO g(8), (32) 
z Jo 

where g(6) is the angular distribution (6 being the angle 
between the outgoing electron momentum and the direc- 
tion normal to the interface) and 

I v = -e 2 VN(0)v F A = J V A , (33) 

with V the applied voltage, A the interface area, and 
vp the Fermi velocity. Eq. I|33|) may be written as 
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ly = (2e 2 /h)NV, where N = Ak F /4w is the number 
of transverse channels that fit in an interface of area 
A. Thus, ly can be interpreted as the current that 
would flow through a transmissive [T(E Z ) = 1 for all E z ] 
normal-normal interface with the same area and subject 



to the same voltage bias. The r 4 dependence of Ins re- 
flects the simultaneous tunneling of two electrons. 

Using the previous definition g z = q z /kp, the angular 
distribution for the current through an arbitrary square 
barrier is 145| 



g{9) = 2cos 3 de 2p ° w[1 - b( - cos9 '> ] 



- [ dg z — 

n Jo \Qz 



cos 2 9) 2 + 5 2 



qI [a(g z , cos0)f e Pow[1 ' b ^ )] 



(34) 



where g c = E c / Ep is the cutoff and the functions a 
and b were defined in Eqs. ((10(1 and 111(1 . 

For S — > 0, we have 



lim 

5^0 



S IT 
d e (g 2_ x 2)2 +(5 2 g/(g) = 2 (35) 



if, as is the case, g c > 1 > |at|. Therefore, in the limit 
5^0, Eq. yields 

g{6) = 2 e 4po^[i-b(c OS e)] [ & ( cos6 ))]4 cos 5q ( 36 ) 
For large barriers (u 3> 1) and finite S we find (with 



oo) 



9(0) 



cos 5 9 + cos 3 6>v/cos 4 6» + <5 2 (37) 
2 cos 5 6» [l + 0(£ 2 )] . (38) 



Combining Eqs. ((32(1 and 1(37(1 . we obtain for the total 
current 



'NS 



1 

12 
1 



IVT 4 



I + (1 + s 2 f/ 2 - s 3 



IyT 4 [1 + 0{8 2 



(39) 
(40) 



However, if the cutoff g c remains finite, Eq. 1(40(1 must 
be replaced by 



T l T ^ 
J-ns — 7^y T 

6 



07T g c 



(41) 



i.e. a finite cutoff qualitatively affects the leading low-S 
dependence of Ins- 

The underlying physics goes as follows. The product 
of hopping matrix elements appearing in 1(31(1 satisfies 

T k * iq T k * 2j _ q oc ku k 2z q 2 S(k n + kaiiWkm - q„) . (42) 

Thus, when crossing the barrier, electrons forming a 
Cooper pair of momenta (q, — q) undergo the following 
process: Their opposite interface-parallel momenta are 
conserved (k^y = qn and k 2 || = — qn). By contrast, 
one of their perpendicular momentum components (more 



specifically, the negative one pointing away from the in- 
terface) is reversed so that both electrons enter the nor- 
mal metal with perpendicular momenta k\ z , k 2z > 0. In 
the limit of eV — > the modulus difference between ki z 
and k\ z is negligible. This means that the electron cur- 
rent through a broad interface will propagate into the 
normal lead in the form of two rays which are symmet- 
ric with respect to the direction normal to the interface. 
Due to axial symmetry, g is only a function of the zenithal 
angle 9 G [0,tt/2]. 

The normalized angular distributions for several bar- 
rier heights are depicted in Fig. [3] in the limit 6 — > 0. 
The lowest barrier which we have considered has u = 1.1. 
This means that, for a typical value of Ep — 5 eV, the dif- 
ference between the height of the barrier and the Fermi 
energy is 0.5 eV, i.e. large enough to ensure that the 
junction operates in the tunneling regime. In Fig. 
finite-height barriers are taken to have a width w = h\p. 
For large Uq we reproduce the analytical cos 5 9 behav- 
ior given in Eq. ((38(1 . As the barrier height decreases, 
the angular distribution becomes more focussed in the 
forward direction because transmission is more sensitive 
to the perpendicular energy. Thus the relative fraction 
of Fermi surface electrons crossing the interface with E z 
close to the highest value Ep increases. That majority 
of transmitted electrons have low parallel momenta and, 
accordingly, a characteristic parallel wave length much 
larger than A p . We will see later that this perpendicular- 
energy selection bears consequences on the length scale 
characterizing the dependence of the total current on the 
radius of the interface. 

In general, knowledge of the current angular distribu- 
tion is physically relevant, as one is ultimately interested 
in directionally separating the pair of entangled electron 
beams for eventual quantum information processing. To 
acquire a more complete picture, we may compare the 
previous results with the case of a NN interface. In that 
case the total tunnel current is 

C7T/2 

(43) 



(44) 



Jnn = Ivt z / d9 sm6 g{6) , 
J o 

where ly is given in Eq. 1-i-SI) and, for large u, 
g{9) = 2cos 3 0. 



9 



Thus we see that electron transport through a tunneling 
NN interface also exhibits focussing which is however less 
sharp than in the NS case [see Eq. Ij38|l ]. The term r 2 in 
Eq. I)44|) reflects the dominance of single-electron tunnel- 
ing at the NN interface. Finally, we may compare Eqs. 
(|38|1 and (|44|l with the cos 9 distribution law exhibited by 
electron current in the bulk of a disordered wire |4o| . 



B. Universal relation between NN and NS 
tunneling conductances 

In the case of a normal interface with high barrier, the 
total current can be integrated to yield 



h 



NN 



I V T 2 



2e2 
h 



-NV. 



(48) 



A. Connection with the multi-mode picture 

We could have derived the angular distributions in Eqs. 
(E2J), GEl, (EH an d (EJ following the scattering theory 
of conduction in normal |47j and normal-superconducting 
|29j | multichannel wires. For an NN interface we can write 
the relation between conductance and transmission prob- 
abilities at the Fermi energy as 



Gnn = 



2e 2 



N 



m-^ _ 2e 2 ^ 



N 



h 



(45) 



n=l 



where A n (n = 1,2,... TV) are the eigenvalues of the 
transmission matrix tt + at the Fermi energy and T n = 
Sm=i l^nm| 2 are the modal transmission probabilities at 
the same energy |2^, which is what we calculate. The 
exchangeability of T\ A n and ^2 n T n reflects the invari- 
ance of the trace |48j. Now consider the transmission 
probability through a square barrier given in Eq. I|15|) . 
We replace T n -> T(E Z ). For E z /U < 1, we have 
T(E Z ) oc Eg oc cos 2 9. Moreover, the sum over transverse 
modes can be replaced by an integration over the zenithal 
angle, J2 n ~^ cast, x J^ 2 d9 sin 9 cos 9. Altogether, the 
angular distribution follows the cos 3 9 law expressed in 
Eq. 0. 

A similar line of argument can be followed for the An- 
dreev current through a NS interface, whose conductance 
is given by 



2e- 



N 

£ 

n=l 



2A 2 



A, 



(46) 



As noted in Ref. 23], the equivalence invoked in Eq. 
(f4"5jl is no longer applicable in Eq. (|47)|) because of its 
nonlinearity. Nevertheless, in the tunneling limit one has 
A n <C 1 and Gns can thus be approximated as 



Gns 



2e^ 
h 



N 



EA 2 Op 2 " T 2 



ti 2 



(47) 



where the second equality is possible because J2 n A n — 
Tr(tt + ) 2 . Arguing as we did for the NN conductance, it 
follows that g(9) oc cos 5 9, which confirms Eq. I|38() . We 
note here that, in Refs. psl l2flj . the Andreev approxi- 
mation was made whereby all the momenta involved are 
assumed to be equal to kp . In our language, this corre- 
sponds to taking S — > in Eq. (|34l) and thereafter. 

Finally, comparison of Eqs. 145|l and l|47(l also illumi- 
nates the contrast between the factor t 4 /2 in Eq. I|32|) 
and the factor t 2 in Eq. (021 . 



Thus t 2 /2 is the average transmission per channel 
[49l l50| . In one dimension (N = 1) one has Inn = 
(2e 2 /h)VT 2 . Eqs. ®, and gHJ) suggest the uni- 

versal relation 



IvIns GvGns 



(Inn) 2 



(Gnn) 2 



(49) 



where G, = h/V (i = V,NS,NN). Eq. O indicates 
that knowledge of Gns and Gnn may allow us to infer Gy 
and, from (|33f) . the effective area of a tunneling interface. 



C. Comparison with the quasiparticle scattering 
method 

Blonder et al. Q studied transport through a onc- 
dimensional NS interface modelled by a delta-barrier one- 
electron potential [U(z) = HS(z)] by solving for the 
quasiparticle scattering amplitudes. If the dimension- 
less parameter Z = mH/h 2 kF is employed to character- 
ize the scattering strength of the barrier, the tunneling 
limit corresponds to Z > 1, for which they obtained 
j ns = !v/Z 4 assuming k B T < eV < A < E F (i.e. 
a low-transmission regime in which Andreev reflection 
is the only charge-transmitting channel). Later, Kupka 
generalized the work of Ref. [3| to investigate the sensi- 
tivity of Andreev and normal reflection to the thickness 
of the barrier [2(]] and to the presence of a realistic 3D 
geometry [2l|. For the case of a broad interface in the 
tunneling limit he obtained J^s = Iv/QZ 4 . Therefore, 
Kupka found a result identical to Eq. (|40|) (to zeroth 
order in 6) with r replaced by 1/Z. In fact, it is easy to 
see that, in the case of a delta-barrier with Z> 1, the 
transparency defined in section IV is precisely r = 1/Z. 
Therefore, comparison of Eqs. (|38fl and 140|) with the re- 
sults of Ref. |2l| completes the discussion of section II by 
establishing the quantitative equivalence between the pic- 
tures of quasiparticle Andreev reflection and two-electron 
(or two-hole) emission. We note that, in Refs. [3I I20I l2l|| . 
the Andreev approximation (5 — > 0) was made. 



VII. CURRENT THROUGH A CIRCULAR 
INTERFACE OF ARBITRARY RADIUS 

In this section we investigate transport through a cir- 
cular NS tunneling interface of arbitrary radius. The 
setup is as depicted in Fig. |2Ja). To make the discus- 
sion more fluent, lengthy mathematical expressions have 
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FIG. 4: Radial dependence of the normalized NS Andreev 
tunnel current through a circular interface of radius R for 
different barrier heights. J(R) = I(R)/tvR 2 and Joo is the 
current density in the thermodynamic limit. Finite barriers 
have a width w — 5Xf- Everywhere A/Ef <C 1 is taken. 
Dots correspond to numerically exact results. Solid lines are 
computed with an approximation described in Appendix B 
which becomes exact for R/Xf 2S> 1. The inset magnifies 
results for small R. 



been transferred to Appendix B, leaving here the presen- 
tation of the main results, which include some analytical 
expressions for the limit of small gap and high barrier. 



Total current 



a factor R 2 to the amplitude, regardless of the incident 
angle. The Cooper pair amplitude becomes ~ R 4 , which 
leads to the R 8 behavior for the probability. 

It is interesting to compare the R 8 law here derived 
with, e.g. the R 4 behavior of the NN tunnel current 
(u^$> 1), namely, 



'-J v r 2 k 2 F R 4 , 



(51) 



or with the i? 6 dependence for the transmission of pho- 
tons through a circular aperture |5ll |. 

Eqs. H50f) and Eq. I|51(l yield the following relation for 
the narrow interface conductances: 



h 
4^2' 



Gns — — Gnn (R 



0) 



(52) 



It is important to note that Eq. (|52|) still applies if 
both conductances are replaced by their momentum- 
independent counterparts. 

In Fig. 01 we plot the current density as a function of 
the interface radius. Dots represent the exact calculation 
taken from Eqs. (|B3(1 and l|B5(l . which we have been 
able to evaluate numerically for u — > oo (up to R — 
1.65Aj?) and u = 10 (up to R = 3X F ), while solid lines are 
obtained from a large-radius approximation described in 
Appendix B. For u = 1.1 convergence problems prevent 
us from presenting numerically exact results. We find 
that the small-radius approximation (~ R 8 ) is correct 



within 1% accuracy up to R 
it overestimates the current. 



O.IAf- Above that value 



The most general expression for the current is given in 
Eq. (jBTJ. Below we focus on the limit S, u 1 < 1. We 
find three regimes of interest, depending on the value of 
R/X F . 

1. Small radius (R <IC Xf) 

This limit is not physically realizable, at least with cur- 
rent materials. However, it is interesting for two reasons. 
First, it yields a radius dependence that directly reflects 
the entangled nature of the electron current. Second, it 
can be used as a unit of current such that, when referred 
to it, calculated currents have a range of validity that 
goes well beyond the geometrical model here considered. 
That permits a direct comparison between different the- 
oretical models and experimental setups. 

For fcpi? < 1 we obtain 

I(R) ~ J = ^J v r 4 k 6 F R 8 . (50) 

This R 8 behavior is easy to understand. To compute the 
current we must square the matrix element between the 
initial and the final state, i.e. the Cooper pair hopping 
amplitude. The tunneling of each electron involves an in- 
tegral over the interface, which for fcpi? <C 1 contributes 



2. Intermediate radius (Xf < R < oo) 

In this region no analytical expression for the current 
is possible. Above R sa 2Xf even the numerical calcula- 
tion of Eq. (|B5|I (which presumes S, u^ 1 <C 1) is difficult, 
since for large radii we cannot compute five strongly os- 
cillating nested integrals. A set of two approximations 
which reduces the number of nested integrals from five 
to three is discussed in Appendix B and expressed in Eqs. 
lHHI and (|B8|) . 

In Fig. 01 we plot I(R)/I(R — > oo), which is the total 
current normalized to the thermodynamic limit expres- 
sion with A in Eq. replaced by irR 2 . For finite 
barriers, w = 5Xf has been taken. A free parameter has 
been adjusted to fit the numerically exact result in the 
region where it is available. As explained in Appendix 
B, such a scheme is particulary well suited for moderate- 
to-large radius values. The inset of Fig. Q] shows that, as 
expected, the approximation fails for small values of i?, 
where it yields an R 4 behavior instead of the correct R 8 
law, thus overestimating the current. 

Here we wish to remark that, unlike in the case of a 
clean NS point contact [23|,|2!|, the radial current depen- 
dence shows no structure of steps and plateaus as more 
channels fit within the area of the interface. This is due 
the fact that we operate in the tunneling regime, which 
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decreases the height of the possible steps and, more im- 
portantly, to the strongly non-adiabatic features of the 
structure along the ^-direction. 

3. Large radius (R — » oo) 

While a numerically exact calculation is already non- 
feasible for R a few times Xp, the approximation de- 
scribed in Appendix B becomes increasingly accurate for 
large R. This allows us to conveniently investigate how 
the broad interface limit is recovered [see Eqs. (|32fl and 
Q40[l] . Such a limit is characterized by I(R) growing with 
R 2 , i.e. proportionally to the area, a behavior also shown 
by the NN conductance. Convergence to the thermody- 
namic limit is much slower for low barriers than for large 
barriers. The reason has to do with focussing. The wave 
length of the characteristic energies E\\ = E — E z deter- 
mines the length scale over which the relative phase be- 
tween distant hopping events varies appreciably. This is 
the distance over which multiple hopping points (which 
play the role of multiple "Feynman paths" ) cancel de- 
structively for large radius interfaces. As discussed in 
the previous section, low barriers are more energy selec- 
tive, making most of the electrons leave with E z close to 
Ep and thus with small E\\ . As a consequence, saturation 
to the large radius limit is achieved on the scale of many 
times Xp . By contrast, high barriers are less energy selec- 
tive and give a greater relative weight to electrons with 
low E z and high En . A large fraction of the electrons 
has a short parallel wave length. This explains why, for 
high barriers, the large R behavior is reached on a short 
length scale. 



B. Length scales in the thermodynamic limit 

It is known that pairing correlations between electrons 
decay exponentially on the scale of the coherence length 
£o = Xp/ir 2 6. This fact is reflected by the exponential 
factors contained in the integrands of the equations for 
I(R) in Appendix B. Thus one might expect that the 
thermodynamic limit relies on such a decay of correla- 
tions. 

The following argument might seem natural. The dou- 
ble integral over the interface of area A may be viewed 
as an integral of the two-electron center of mass, which 
yields a factor A, and an integration over the relative co- 
ordinate, which is independent of A due to a convergence 
factor which expresses the loss of pairing correlations. 
The final current would grow as / ~ A£q ~ A/5 2 . How- 
ever, as discussed in the previous subsection, the ther- 
modynamic limit is achieved on a much shorter scale, 
namely, the Fermi wave length. If an electron leaves 
through point i"i one may wonder what is the contri- 
bution to the amplitude stemming from the possibility 
that the second electron leaves through r2, eventually in- 
tegrating over Y2- Eq. 124|) suggests that the amplitude 




FIG. 5: Angular correlation profile (in arbitrary units) of the 
conditional probability distribution P(fi|fio) that, in a given 
tunneling event, an electrons goes into Q if the other electron 
has gone into Qo- Here we plot P(f2|fio) as a function of £1 
for fixed £7o = (0o,4>o) = (7r/4, 0). From top-left to bottom- 
right the radii are: R — 3, 1, 0.5, O.OIAf- Observe that, as 
R increases, the angular dependence of the second electron 
tends to be the conjugate of the first one, i.e. the distribution 
becomes peaked around Q — (7r/4, n). Note also that, for 
small R, P(Q\Qo) becomes oc cos 2 9 regardless of Qo- 

for two electrons leaving through ri and r2 will involve 
the sum of many oscillating terms with different wave 
lengths, the shortest ones being ~ Xp. This reflects the 
interference among the many possible momenta that may 
be involved in the hopping process. Such an interference 
leads to an oscillating amplitude which decays fast on the 
scale of Xp, rendering the exponentially convergent factor 
irrelevant. Thus, in the thermodynamic limit the current 
tends to a well defined value for ^ -> oo (i5 -> 0). In Ap- 
pendix B we provide a more mathematical discussion of 
this result. 

One may also investigate the first correction for small, 
finite 5. As indicated in Eqs. I|39|) and (|4(jp. it increases 
the current. However, in the presence of a finite cutoff 
(g c < oo), a nonzero value of 5 generates the opposite 
trend. As discussed in Appendix B, at tiny relative dis- 
tances between hopping points (|ri — ra| < 5Xp), the 
amplitude increases considerably. A finite upper momen- 
tum cutoff rounds the physics at short length scales, thus 
eliminating such a short-distance increase. The result is 
that, with a finite cutoff, the first correction to the 5 = 
limit is a decreasing linear term in 5, as revealed in Eq. 

C. Angular distribution and correlation 

We have computed the conditional probability dis- 
tribution P(f2|f2o) for an electron to be emitted into 



12 



il = (9, 4>) given that the other electron is emitted in 
a fixed direction Q . Such a distribution is shown in 
Fig. El for £!o = (7r/4,0). We observe that, for large 
R/Xf, the angular distribution of the second electron 
is quite focussed around Q = (tt/4, 7r), which is mirror- 
symmetric to Hq. As R/Xf decreases, the angular cor- 
relation between electrons disappears and, as a function 
of il, P(f2|f2o) becomes independent of the given value of 
In particular it tends to ~ cos 2 9. 
We may also study the probability distribution that 
one electron is emitted into direction Q regardless of the 
direction chosen by the other electron. This amounts 
to the calculation of an effective g(9) for a finite radius 
interface to be introduced in an equation like Eq. I|32[l to 
compute the current (by symmetry, such a distribution is 
independent of cj>). As expected, one finds such effective 
angular distribution to be ~ cos 5 9 for large R [see Eq. 
ijSSJ], which contrasts with the sharp ^-dependence of 
the conditional angular distribution P(f2|Qo) for given 

n . 

For small R, the effective g{9) goes like cos 2 9, i.e. it 
becomes identical to P(f2|f2o)- This coincidence reflects 
the loss of angular correlations. The cos 2 9 behavior may 
be understood physically as stemming from a random 
choice of final k|| , which yields a cos 9 factor (since |k|| | = 
fcp sin#), weighted by a cos 9 reduction accounting for 
the projection of the current over the z direction. An 
equivalent study for a NN interface yields also g{9) oc 
cos 2 9. Thus we see that the loss of angular correlations 
after transmission through a tiny hole makes the NN and 
NS interfaces display similar angular distributions. 

The crossover from g{9) oc cos 2 9 to cos 5 9 as R in- 
creases involves a decrease of the width A9 of the angu- 
lar distribution. A detailed numerical analysis confirms 
this result but reveals that A9 is not a monotonically 
decreasing function of R (not shown). 



VIII. NONLOCAL ENTANGLEMENT IN A 
TWO-POINT INTERFACE 

Let us turn our attention to a tunneling interface con- 
sisting of two small holes, as depicted in Fig. Hfb). By 
"small" we mean satisfying R/Xp <C 1 . This is the limit 
in which the detailed structure of a given hole is not im- 
portant and the joint behavior of the two holes is a sole 
function of their relative distance r and the current that 
would flow through one of the holes if it were isolated. 
We expect the conclusions obtained in this section to be 
applicable to similar interfaces made of pairs of different 
point-like apertures such as, e.g. two point-contacts or 
two quantum dots weakly coupled to both electrodes § . 

The current through a two-point interface has three 
contributions. One of them is the sum of the currents 
that would flow through each hole in the absence of the 
other one. Since the two orifices are assumed to be iden- 
tical we refer to it as 2Iq, where 1$ is given in Eq. 1)50(1. 
This contribution collects the events in which the two 
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FIG. 6: Current density through a two-point interface stem- 
ming from non-locally entangled electron pairs, as a function 
of the distance between points. The dashed line corresponds 
to the current obtained using an energy-independent hopping 
approximation whereas in the solid line the correct momen- 
tum dependence has been taken into account. 



electrons tunnel through the same opening. A second 
contribution I e (r) comes from those events in which each 
electron leaves through a different hole. This is the most 
interesting contribution since it involves two non-locally 
entangled electrons forming a spin singlet. The third con- 
tribution, Iiir), accounts for the interference between the 
previous processes. 
If we write 



I = 2h + I e {r) + h{r) 



(53) 



we obtain for the entangled current in the high barrier 
limit 

J B (r) - l8I [B 2 (k F r) + F 2 {k F r)B 2 (k F r)] , (54) 
where B(x) in Eq. (|B10|) and 

sin x — x cos x 



F{x) =3 



(55) 



For 5 <C 1, and noting that we are not interested in tiny 
distances r < SXf, we can write 

I e (r) = 2I Q [F 2 (k F r) + F 4 (k F r)}e- 2r/7rio . (56) 

This is a fast decay because of the geometrical pref- 
actor, which goes like r for k F r 3> 1. For instance, 
I e (£ o )/I e (0) ~ 10" 15 , with data taken from Al (£ ^ 
10 3 Af)- For possible comparison with other tunneling 
models it is interesting to write the entangled conduc- 
tance G e (r) = I e (r)/V in terms of the normal conduc- 
tance through one narrow hole, Gnn- Using Eq. l(52"|l . 
we obtain 



Ge(r) = ^G 2 ™ [F 2 (k F r) + P 4 (^r)] e - 2r /^ . (57) 
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To keep track of the interference terms, it is convenient 
to adopt a schematic notation whereby = t a + t& is 
the tunneling Hamiltonian through points a and b. Then 
one notes that, as obtained from Eqs. ((25(1 . ((26(1 and 
(|29|l . the total current can be symbolically written as I ~ 
\(ta + t b )(t a +t b )\ 2 . In this language J ~ \t a t a \ 2 - The F 2 
term in l(56J) corresponds to ~ |£ a £b| 2 + |Wa| 2 , while the 
F A term stems from the interference ~ (t a tb)(tbt a )* +c.c. 
Altogether, I e (r) ~ \t a t b + t b t a \ 2 . 

The interference current may be divided into two con- 
tributions, 



Ii = hi 



HI i 



(58) 



corresponding to the different types of outgoing chan- 
nel pairs which may interfere. The first contribution 
stems from the interference between both electrons leav- 
ing through point a and both electrons leaving through 
point b, In ~ {t a t a )(tbtb)* + c.c. One obtains 



I il (r) = 2I F 2 (k F r) 



(59) 



Hi{t) comes from the interference between the channel 
in which the two electrons leave through the same hole 
and that in which they exit through different openings, 
-f»2(r) ~ (t a ta)(t a tb)* + c.c, plus three other equivalent 
contributions, altogether summing 



Ia{r) = 8I F 2 (k F r)e~ r ^° 



(60) 



In the hypothetical case where orifices a and b are 
connected to different normal electrodes [e.g. when an 
opaque barrier divides into two halves the normal metal 
of Fig. Etb)], the interference contributions l(59() and 160(1 
would be absent. Then one would have I = 2Iq + I e (r). 



IX. FAILURE OF THE 
MOMENTUM-INDEPENDENT HOPPING 
APPROXIMATION 

It has been common in the literature on the tunneling 
Hamiltonian to assume that the tunneling matrix ele- 
ments appearing in Q are independent of the perpen- 
dicular momenta k z q z (see, for instance, Ref. |35j). Be- 
low we show that, for three-dimensional problems, such 



an assumption is unjustified and leads to a number of 
physical inconsistencies [52| . 

For simplicity we focus on the high barrier limit. 
To investigate the consequences of the momentum- 
independent hopping approximation, we replace Eq. 1(121) 
by 



kq 



nN(o) 



<5(k|| -q.\\)k 2 F 



(61) 



i.e. we change k z q z by k F . 

Broad interface. For a large NS junction, we find that 
the total current in units of Iy diverges (x = cos 9): 



h 



NS 



r 4 Iy 



1 dx x 2 + \Ar 4 + S 2 



oo . 



(62) 



x 2(x 4 + S 2 ) 
i.e. Ins grows faster than A for A — > oo. Eq. ((62(1 is the 
analogue of Eqs. l(52j) and (|3T|) . 

A different divergence occurs for a broad NN tunnel 
junction: 



Inn = 2r ly 



1 dx 



oo , 



(63) 



J^dxx 3 



which contrasts with the finite integral 7nn 
obtained from inserting 1(44(1 into 1431) . 

Local Hamiltonian. If one attempts to derive the real 
space tunneling Hamiltonian with the assumption 161II . 
one obtains an expression identical to that in Eq. I(19() 
with L(z, z') replaced by 



M(z,z') 



(64) 



' k z ,q z 



As in section IV, we use stationary waves for ipk z .q z (z) 
Invoking the identity 



Esin(k z z) =V — 



(65) 



k z >0 



we obtain 



J 



2n 4 N(0) 



f , f° dz f L dz' , t . . , . , . 
dr — —ip , N {r,z;a)^ s {r,z;a)+R, 

J J-L z Jo z 



(66) 



r 



where the reference to the principal value has been re- ((64(1 . we would have obtained a different Hamiltonian, 
moved because, in the tunneling limit, the fields vanish 
linearly at the origin. 

If we had chosen plane wave functions for fk z ,q z hi Eq. 
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namely, 



^ 8ir 2 N(0) 



dr^ N (r,0;a)ip s (r,0;a)+R.c. 



(67) 

which is some times proposed in the literature (see e.g. 
Ref. 0). This situation, whereby plane- wave and 
stationary- wave representations lead to different, both 
unphysical, local Hamiltonians contrasts with the sce- 
nario obtained with the right matrix element. As noted 
in section IV, the more physical choice (|12f) leads in both 
representations (plane and stationary waves) to the cor- 
rect local Hamiltonian The fact that Eq. l|6*T|) 
leads to a wrong real space Hamiltonian which, more- 
over, depends on the choice of representation, may be 
viewed as further proof of the inadequacy of the energy- 
independent hopping model. 

Thermodynamic limit. For a NS interface with d — > 0, 
a dimensional analysis for A — ► oo suggests that the total 
current 7ns diverges non-thermodynamically like ~ A 2 . 
For a NN interface, we find the divergence A In A. 

Unitarity. The divergences expressed in Eqs. (|62H and 
(|63|l . as well as the related anomalous thermodynamic 
behavior, could have been anticipated by noting that, if 
Tk q is assumed to be independent of energy, then Eq. 
(TH)|l must be multiplied by (E F /E Z ) 2 . As a result, the 
transmission probability at energy E z , which should stay 
smaller than unity, grows instead as T(E Z ) ~ E^ 1 for 
E z — > 0. Such a violation of unitarity necessarily gener- 
ates a divergent current in the broad interface limit for 
both NN and NS interfaces. 

Nonlocally entangled current. Finally, we note that, 
using (|rTT]) , the nonlocally entangled current through two 
distant points is 



I e {r) =2I [F 2 (k F r)+F 4 (k F r)]e~ 2r /^ , 



where 



F(x) = 



SIM 



(68) 



(69) 



with the tildes generally referring to the momentum- 
independent approximation. Here, Iq = 81Iq is the cur- 
rent through one narrow hole. Correspondingly, the en- 
tangled conductance G e (r) is written like in Eq. (|57|> 
with F(kpr) replaced by F(k F r). 

Comparison of Eqs. I|56|l and l|68|l indicates that the 
r-dependence of the geometrical prefactor is markedly 
different: For growing r, the nonlocally entangled cur- 
rent decays much more slowly (r~ 2 ) than its momentum- 
dependent counterpart {r~ 4 ). It is interesting to compare 
the ratios A(r) = 7 e (r)/J e (Q) and A(r) = / e (r)/I e (0). 
While A(0) — A(0) = 1 by construction, the ratio A/A 
becomes ~6x 10~ 4 and 2 x 1CT 7 for r/\ F = 20 and 10 3 , 
respectively. 

Interference terms. As expected from the comparison 
of Eqs. (|56|l and l|68|l . the interference contributions are 
identical to those discussed in the previous section with 
F(k F r) replaced by F(k F r) in Eqs. (£>HJ) and (JHOJ). 



Generality of the model. An important question is 
whether our results for the entangled and interference 
current through pairs of tiny geometrical holes apply to 
other, more realistic pairs of small interfaces such as two 
point contacts or two quantum dots |g. The fact that 
the decay with distance of the entangled current reported 
in Refs. 0, 0, IE E3, EE EE E3, Ell follows the same 
law as Eqs. (|68l) and (|69|) (except for the F A term there 
neglected), suggests that such is indeed the case. Below 
we prove this expectation. 

Due to Eq. the sum in Eq. lj3T|) involves 



E n 



-q z e 



(70) 



This sum over q is clearly affected by the presence of the 
q 2 factor, yielding a result cx F(k F r), with r = |ri — r 2 |. 
In the momentum-independent hopping approximation, 
q 2 is replaced by k F , rendering the sum cx F{k F r). In 
fact, the two functions are related: 



dz 2 



F(k F \Jr 2 + z 2 ) 



= \k 2 F F{k F r) . (71) 



We note that the distance dependence is determined by 
the properties of the superconductor and not by those 
of the normal electrode. If a quantum dot mediates be- 
tween the superconductor and the normal metal, then 
an effective hopping must be introduced in (|70|l which, 
however, does not add any new momentum dependence 
[see Eq. (11) of Ref. Departure from the specific 

type of contact here considered will translate only into 
a different value of Joj the distance dependent prefactor 
remaining identical. We notice, however, that the preced- 
ing discussion is restricted to the case where quasiparticle 
propagation is ballistic in both electrodes, i.e. we neglect 
the effect of impurities or additional barriers [54| . 



X. SUMMARY 

We have investigated the electron current through a NS 
tunneling structure in the regime fc^T < < A where 
Andreev reflection is the dominant transmissive channel. 
We have rigorously established the physical equivalence 
between Cooper pair emission and Andreev reflection of 
an incident hole. A local tunneling Hamiltonian has been 
derived by properly truncating that of an infinite inter- 
face in order to describe tunneling through an arbitrarily 
shaped interface. Such a scheme has been applied to 
study transport through a circular interface of arbitrary 
radius and through an interface made of two tiny holes. 
In the former case, the angular correlations between the 
two emitted electrons have been elucidated and shown to 
be lost as the interface radius becomes small. We have 
also investigated how the thermodynamic limit is recov- 
ered, showing that, due to the destructive interference 
between possible exit points, it is achieved for radii a few 
times the Fermi wave length. For the case of a two-point 



15 



interface, we have calculated the nonlocally entangled 
current stemming from processes in which each electron 
leaves the superconductor through a different orifice. We 
have found that, as a function of the distance between 
openings, such an entangled current decays quickly on 
the scale of one Fermi wave length. The interference be- 
tween the various outgoing two-electron channels has also 
been investigated and shown to yield contributions com- 
parable to the nonlocally entangled current. We have 
found that, in a three-dimensional problem, it is impor- 
tant to employ hopping matrix elements with the right 
momentum dependence in order to obtain sound physical 
results in questions having to do with the local tunneling 
Hamiltonian (whose correct form has also been obtained 
from a tight-binding description), the thermodynamic 
limit, the preservation of unitarity, and the distance de- 
pendence of the nonlocally entangled current through a 
two-point interface. An important virtue of the method 
here developed is that it enables the systematic study of 
Cooper pair emission through arbitrary NS tunneling in- 
terfaces and opens the door to a convenient exploration 
of the fate of Cooper pairs in the normal metal and, in 
particular, to the loss of phase and spin coherence be- 
tween emitted electrons. 
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APPENDIX A: DISCRETE VS. CONTINUUM 
SPACE 

Take a discrete chain made of N sites with period a 
described by the Hamiltonian 



N-l 



(Al) 



i=i 



where t = h 2 /2ma 2 > is the hopping parameter that 
yields an effective mass m in the continuum limit. 
The eigenstates of this chain are of the form 



N - 



1/2 N 



- j ^2 sin(fc„2i)cj |vac) , (A2) 



where z, = ia and k n 
The eigenvalues are 



im/a(N + 1) with i,n6 [1, N]. 



E n = -2tcos(fc„a) . (A3) 
The basis set {|</> n )} is orthonormal. Thus we may write 

1/2 N 




-j ^sin(fc„^)cj» (A4) 



i=l 
1/2 N 



-J ^sin(fc„z. i )4„ • ( A5 ) 



We write the transfer Hamiltonian between two A-site 
chains as 



-t'a[b^x + H.c 
2f N N 



N ■ 



s ' m (k n a) sm(k m a)al n b k 



(A6) 
H.c.(A7) 



which may be treated as a small perturbation when t' <C 
t. 

To investigate the continuum limit, we take a — ► and 
t — > oo so that m and kp remain finite. We also take 
N ^> 1. Noting that the sine functions in l|A7(l can be 
approximated by their arguments k n a < kpa <C 1 and 



that k„ 



iv /a — > co, we get 



Ht — 



2t'a 3 



L 



^ kqa\b q + H.c. 



(A8) 



fc,g>0 



This Hamiltonian is bilinear in the momenta of the elec- 
tron on the right and left chain. If we were in 3D we 
would specify that the bilinearity refers to the momenta 
perpendicular to the interface plane. This hamiltonian is 
analogous to that which we proposed for the continuum 
Bardeen theory in the case of a high barrier [see Eqs. (@J 
and H2J)]. 

We may work out the corresponding Hamiltonian in 
real space. For that we note that, in the continuum limit, 
Ht in Eq. (|A6|I can be expressed in terms of field op- 
erators evaluated in z = ±o. When a — > 0, the field 
operators can be expanded as 



dipiz) 



dz 



+ 0(a 2 ) , (A9) 



2 = 



where ip(0) = is a condition that results naturally from 
the properties of the wave functions in a chain starting in 
i = 1 or i = — 1. For such chains, i = is an imaginary 
point where the wave function necessarily vanishes; it is 
the place where we would locate the hard wall in a con- 
tinuum description |55j . Then the tunneling Hamiltonian 
can be written 



Hn 



t l a 3 d ^ R ( Z ) 



dz 



d^ L {z') 



dz' 



H.c. 



(A10) 



16 



This Hamiltonian is exactly the one-dimensional version 
of that in Eq. \2'.\\) . The fact that we have derived it from 
a completely different set of physical arguments should 
be viewed as a definite proof of the adequacy of the tun- 
neling Hamiltonians proposed in section IV. The Hamil- 
tonians IjASjl and (|A10|I have been obtained in the con- 
tinuum limit. On the other hand, Eqs. (|12|) and l|23l) 
were derived for high barriers or, equivalently, low ener- 
gies. Clearly, this is not a coincidence, since it is at low 
energies where the long wavelengths make the electron 
move in the chain as in continuum space. 

APPENDIX B: TOTAL CURRENT VS. 
INTERFACE RADIUS 

To calculate the total current as a function of the in- 
terface radius R we have to evaluate the matrix element 
(|31|l using hopping energies obtained from the tunneling 
Hamitonian (|24[1 . In the resulting expression we need 
to integrate over the final momenta of the two electrons 
in the normal metal, the momentum of the intermediate 
virtual state consisting of a quasiparticle in the super- 



conductor, as well as the coordinates of the points where 
each electron crosses the interface area. The integrations 
over the momenta in the final state lead to four angular 
integrals ($1,2 S [0, n/2]; ipi : 2 & [0,27r]), the moduli being 
fixed by the condition fcsT, eV — > 0. The integration over 
the superconductor excited states leads to three integrals: 
9 S £ [0, n],ip s 6 [0, 2n], q S [0, q c \. On the other hand, in- 
tegration over the hopping points of each electron leads 
to two interface integrals (n,2 6 [0, R], Q r \. r 2 € [0, 27r]), 
which makes four more integrals, totalling eleven real 
variables to be integrated. Using the symmetry prop- 
erty that the integrand is independent of one azimuthal 
angle, and solving analytically the four real space inte- 
grals, we are left with six non-reducible nested integrals 
of strongly oscillating functions. 

We define k = k/fc^?, g = q/fci?. Since the mod- 
ula of the final momenta are fixed by conservation re- 
quirements, we may write K\u — sin#;, n Z i — cos#; 
(i = 1,2). For the virtual states in the superconductor: 
g\\ = £>sin6> s , g z = gcos9 s . 

The general, exact formula for the total current as a 
function or R is 



I(R) 



4 (k F R) 2 f f 2 



2 „2 P0 U'[l-b(Kl i ,K2 i )] 

2z e 



dg 



W *{g 2 -l) 2 + 5 2t 



f dn s e^- b ^ n o(/ S , is fo)j(|^|-/e 1 |,|,|§ii + /s a || 



-, 2 



J = l,2 



where J is a short-hand notation for 



_ J 1 (k F Rx)J 1 (k F Ry) 

J ( x j y) — 

xy 



(Bl) 
(B2) 



The first-order Bessel functions result from the exact integration over the tunneling points ri and r2. 

For 5^0, the Lorentzian becomes a delta function and the integral over g is evaluated exactly. We get (with u 
still arbitrary) 



I(R) = I V T 



(kpRy 

47T 3 



dfti / d K2z K 2 lz K 2 2z e 2pow[1 - a{Kl *' K2 * )] 



fdn s e^- b ^ J] o(^, e ,)J(|^||-^|,|,|^|+^|||) 

« a in 



For S arbitrary and u ^> 1, Eq. (|B1(I becomes 



dQi / dhl2 z Ki z K 



2 „2 
2z 



3=1,2 



dg 



W {g 2 - l) 2 + 6 2t 



g A / dQ 3 J(\g\\ - %x\\\, \g\\ + k 2 \\\) 



Finally, for both 6 — > and 11 > 1, we obtain 

\2 



i{ R ) = i V T i{ ^- J dn,J dK 2z K 2 lz K 2 2z J dn s j(\ 



Q\\ - + «2|||) 



(B3) 



(B4) 



(B5) 



which for kpR <C 1 leads to Eq. 



(|5L)fl in the main text. This is easy to see considering that 
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lim^^o J\{kpRx) / x = kpR/2. 

Even after making <5, u^ 1 —> 0, the resulting expression 
(IB5|) is such that a numerical integration for arbitrary R 
is not yet possible. In order to evaluate (|B3|) and JB5I) 
numerically we need to introduce a set of two approxi- 
mations which are good for k F R 1 and reasonable for 
intermediate R. To introduce the first approximation we 
go back to the original expression where the space 
coordinates have not yet been integrated. Then we shift 
from the two space coordinates (ri,^) to center-of-mass 
and relative coordinates (r c ,r). The integration domain 
of the center-of-mass coordinate r c is still a circle of ra- 
dius R. However, the integration region of the relative 
coordinate r is more complicated: It is eye-shaped and 
centered around r c . The first approximation consists in 
assuming that, for all r c , the integration domain of the 
relative coordinate is circular instead of eye-shaped. The 
area of such a circular region is a free parameter which 
can be adjusted by, e.g. comparing the approximate re- 
sult with the exact calculation for those values of R for 
which the latter can be performed. 

It is intuitive (and rigorously proved in subsection 
VII. B) that, because of diffraction, when R < A^, the 
parallel momentum is not conserved and, in particular, 



the two electrons do not leave necessarily with opposite 
parallel momenta [see Fig. ©]. Nevertheless, as R in- 
creases the interface begins to be large enough so as to 
permit parallel momentum to become better conserved. 
A quasi-delta function <5(kx || -J— k 2 1| ) effectively appears. In 
particular we have: lim^^oo Ji(fc||i?)/fc|| = 27r<5(k||)/i?. 
Thus, our second approximation consists in assuming 
that, for all R > X F , the quasi-delta is an exact delta: 
8 — > 5. This is equivalent to the assumption that there is 
no diffraction, i.e. that we work in the ray optics limit. 
This approximation becomes exact as R — > oo and it is 
a reasonable one for finite radii. Of course, this approxi- 
mation fails for R < Af, yielding a wrong R 4 behavior. 

With the two previous approximations we can reduce 
the number of numerical integrals from five to three. To 
write the resulting expressions, let us introduce some 
compact notation. We define x = cos 9 (where 9 is the 
angle formed by the outgoing momentum with the di- 
rection normal to the interface) , y = cos 9 q (9 q having a 
similar definition within the superconductor), A = fci?|r c |, 
and /j, = k F \r\. 

For 6 — > and arbitrary u we obtain 



I{R)~I v t a / dxx^e 2 ^ 1 -^ 



k F R 2A 



dy- 



2 2 Po w[l-b(y)] 



(k F R) 2 J x 2 -y 21 

VI - 7 j Q (r^X)y/l-x^ .h(r(\)^l^) - ^Jl-x 2 J (r(A)v^7) J^r(A) V~l - z 2 )] }' 



(B6) 



where r(A) is the radius of the approximate circular do- 
main over which the relative coordinate r is integrated. 
If the circle is assumed to have the same area as the eye, 
we obtain 



r(A) 



(k F R) 2 arccos 



A 



k F R 



— A\J (kpR) 2 — A 2 



1/2 



(B7) 

but in practice this criterion is found to overestimate the 
total current. Thus we decide to adopt the ansatz 



r(A) = 2k F R 1 



A 



k F R 



(B8) 



where a is a parameter to be adjusted by comparison 
with the exact solution in those cases where it can be 



computed. In particular, a has been adjusted from 
the last two exact numerical values of each curve, i.e. 
from the two largest computationally possible radii. We 
note that both l|B7|) and (|B8|I satisfy the requirement 
r(A) — > 2k F R for A — * 0. The value a — 1 corresponds 
to the case where the circle is chosen to be the maxi- 
mum circle which fits within the eye-shaped integration 
domain. As expected, this criterion underestimates the 
current. The formula (|B7|I . which overestimates the re- 
sult, can be approximated with a ~ 0.7. Thus it comes 
as no surprise that the value of a obtained by comparing 
with the exact result (when available) is an intermediate 
number, namely, a — 0.84, which has been used for the 
NS curves in Fig. 0] 

For arbitrary S and u ^> 1, the total current becomes 



I(R) ~ I v r 4 



dx x 



(k F R) 2 



k F R 



dXX 



r(A) 



n 2 



d/j, hJq{h^/ 1 — x 2 )B(n) 



(B9) 
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sm{S{5)^)] 



S 2 cos[arctan<5/2 + S(5)fj] 1 -nts)/. 

r 



(BIO) 



r 



where 



S(5) = 



D(5) = 




Thus, for (5 < 1 we may write 

sin(/i) cos(fi + 5/2) 



B(n) 



-Sn/2 



The effect of the phase-shift 5/2 is only appreciable 
for ^ < S, i.e. for r < SXp -C Aj?, as can be seen by 
expanding £>(/i) for small fi: 



The phase-shift generates a divergence for \i — > 0. Al- 
though integrable thanks to the multiplying factor in 
Eq. (|B9J) . this divergence affects the final result. Its 

(BH) 

range of relevance may be estimated by making 5/2fi 
equal to the limiting value 1/3 which one would obtain 
with 5 = 0. This yields a range ro = (3/4tt)5Xf, which 
(B12) will be washed out by any realistic momentum cutoff 
q c ~ kp <C kp/5. 



(B13) Finally, we note that comparison of Eqs. (|56|1 and 

(|B13p clearly reveals that the entangled current I e (r) 
given in (|56|l is essentially proportional to B 2 (kpr). As 
discussed in Sec. IX, I e (r) decays faster than the prefac- 
tor obtained from momentum-independent hopping ma- 
trix elements [see Eq. H68fl ]. The current increase which 
results from such an unphysical approximation translates 
(B14) into a divergent thermodynamic limit (see also Sec. IX). 
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